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ABSTRACT 

A two step hybrid perturbation-Galerkin technique is applied to the problem of deter- 
mining the resonant frequencies of one- or several-degree(s)-of-freedom nonlinear systems 
involving a parameter. In step one, the Lindstedt-Poincare method is used to determine 
perturbation solutions which are formally valid about one or more special values of the 
parameter (e.g. for small or large values of the parameter). In step two, a subset of the 
perturbation coordinate functions determined in step one is used in a Galerkin type approx- 
imation. The technique is illustrated for several one-degree-of-freedom systems, including 
the Duffing and van der Pol oscillators, as well as for the compound pendulum. For all of the 
examples considered, it is shown that the frequencies obtained by the hybrid technique using 
only a few terms from the perturbation solutions are significantly more accurate than the 
perturbation results on which they are based, and they compare very well with frequencies 
obtained by purely numerical methods. 
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1. INTRODUCTION 


In this paper we present and discuss a two-step hybrid perturbation- Galerkin technique 
for the computation of the resonant frequencies of nonlinear oscillating systems with a finite 
number of degrees of freedom involving a scalar parameter e. In previous papers, we have 
developed and applied different versions of the hybrid technique to several classes of two- 
point boundary- value problems for ordinary differential equations [2,6,7], to some boundary 
value problems for elliptic partial differential equations [8], and to some integral equations 
of the first kind which arise in slender body theory [5]. For each of these classes of problems, 
the method has yielded results which are typically far more accurate than the perturbation 
solutions on which they are based, and often provide at least reasonable solutions even for 
values of the expansion parameter for which the perturbation solution is meaningless. 

The idea of exploiting perturbation expansions in conjunction with Galerkin or variational 
techniques was introduced by Noor and Peters in 1979 [12] and developed by Noor and his 
collaborators in a number of papers (see e.g. [11], as well as many other references cited in 
[5]). Noor’s “reduced basis method” is a combination of finite element or other discretization 
techniques, perturbation expansions, and Galerkin (or variational) techniques. 

In general terms, our two-step hybrid technique consists of computing a few terms in the 
perturbation expansion of the solution about one or more values of the parameter e and then 
using a subset of these functions, with new amplitudes, in a Galerkin type approximation. 
This manner of combining the perturbation and Galerkin approaches (which we will describe 
in more detail below) seems to overcome some of the drawbacks associated with each of the 
methods when they are applied by themselves, while preserving some of the good features of 
each one. In particular, the perturbation method has at least two significant drawbacks. The 
first is that, for most practical problems, only a few terms in a perturbation expansion can 
be computed because of the rapidly increasing amount of “labor” that is required to compute 
each additional term. Secondly, the expansion parameter must usually be restricted to values 
which lie close to the point about which the expansion was constructed, in order to obtain 
approximations of acceptable accuracy. A drawback of the Galerkin method is the problem, 
from a practical point of view, of selecting a small number of “good” basis functions. As we 
shall demonstrate below, the functions determined by the perturbation method appear to be 
very effective basis functions and hence our method overcomes the main drawback associated 
with the Galerkin method. Also, the new amplitudes determined by the method produce 
approximations which are typically much more accurate than the perturbation solutions on 
which they are based. 

In the following sections we describe our method in the context of the problem of de- 
termining the resonant frequencies of nonlinear oscillating systems. For simplicity and in 
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order to illustrate explicitly some of the salient features of the method, we consider first 
systems with only one degree of freedom. The method is described in detail for such sys- 
tems and then applied to several model one- degree-of- freedom problems. This allows us to 
obtain several explicit results and expressions, which provide insight into the application of 
the method to more complicated systems. We then describe the method in the context of 
a general (conservative) system with a finite number of degrees of freedom which oscillates 
about a stable equilibrium state. The method is then applied to the classical problem of a 
compound pendulum. Observations about the method are presented in the final section. 


2. ONE-DEGREE-OF-FREEDOM SYSTEMS 


We consider first the problem determining the resonant frequency v of a one-degree-of- 
freedom nonlinear oscillator which we write in nondimensional form as 


V 1 u" + u + (e/ a) f(a u, a v u , e) = 0, 

it(0) = 1, u r (o) - o, 

u(x + 2 it) = u(x) for all x > 0, 

r2 7T 

/ /(a u y a v u , e) u* dx = 0. 

Jo 


(la) 

(lb) 
(1 c) 

( 2 ) 


In Eq (la), a is the maximum amplitude of the response, e is a “small” parameter, and 
/ is a specified nonlinear function of its arguments. (Here the dimensional response y of the 
oscillator has been nondimensionalized by defining u = y/a. Also, the primes in Eqs (l)-(2) 
denote differentiation with respect to the nondimensional time x = tot, where t is time, 
u> = 2 7 r/T is the (unknown) resonant frequency of the oscillator, and T is the (unknown) 
period of the oscillation. Then v = w/w 0 , where w 0 is the linear natural frequency of the 
oscillator.) Conditions (lb) and (lc) insure that u has its maximum amplitude at x = 0 and 
is 2n periodic. Condition (2) is a necessary and sufficient condition for Eq (1) to have a 2i r 
periodic solution. We note that in the special case when f does not depend explicitly on 
u', i.e. / = f(au, e), then Eq (la) represents a conservative system. Then condition (2) is 
satisfied automatically and Eqs (1) can be used to express v as 


«/ = K{f[F(u,e)]- 1/2 du}~\ 

2 e f 1 

F(u,e) = 1 - u 2 + — / f(ax, c) dx, 

Cl Ju 


( 3 ) 


where c is defined by the requirements that c > 0 and F(—c, e) = 0. Thus, for this special 
case, the computation of v reduces to a quadrature, which will be useful for comparison 
purposes with our perturbation and hybrid results (see [9]). 
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3. THE HYBRID PERTURB ATION-GALERKIN TECHNIQUE 


We now apply a slightly modified version of a two step hybrid perturbation-Galerkin 
technique which has been introduced and discussed in a series of papers [2,5-8]. In step one 
of the method, we use the Lindstedt-Poincare method [10] to find approximate solutions for 
u, v, and a in the form 

u = X Uj(x)e 3 + 0(e N ), 

i = o 

" = i + £ »i £> + o( £ "). (4) 

j= 1 

a = X a i eJ + °^ N ), 

3=0 

which are formally valid as t — * 0. In (4), the {uj(x)} are the unknown perturbation co- 
ordinate functions, while the constants {j/j} and {<*_,}, which determine the perturbation 
approximations to the frequency and initial amplitude, respectively, are also unknown. To 
determine these quantities, we substitute the expansions (4) into Eqs (1), expand the result- 
ing equations in power series in e about e = 0, and then equate the coefficients of like powers 
of e on each side of these expressions. In this manner, we obtain a sequence of problems to 
solve for each of the unknowns. In particular, we find easily that 

u 0 (x) = cos(x), (5) 

while each of the functions Uj(x) with j > 1 satisfies a problem of the form 


«" + Uj = 2i/,-cos(x) + gj, 


(6) 


with Uj(0) = u'(0) = 0 and Uj(x + 2w) — Uj(x). Here gj depends on Uk, and a,k with 
k < j. Thus, in order for Uj to be 2x periodic, the right side of (6) must be orthogonal to 
both cos(x) and sin(x), since otherwise “secular” terms proportional to x sin(x) and x cos(x) 
will appear in the solution for itj(x). These conditions yield the relations 


\ r 2ir r 2 7T 

Vj = — - — / gj cos (x)dx, / gj sin(x)dx = 0, 

2 7T JO Jo 


which are two equations for the two unknowns Vj and In particular, for j 

become 

1 f 2n 

;~ 7 ra” Jo ^ a ° COS ^’ ~ a ° s ^ n ( :r )5 0) cos (x)dx, 




2 7r a 0 

r2 7r 


J rl 7 r 

' f(a 0 cos(x), —a 0 sin(x), 0) sin(x)dx. 
0 


( 7 ) 

1 Eqs (7) 

( 8 ) 
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The second equation in (8) is an equation for a 0 , while the first equation expresses Ui as a 
function of do . 

In the examples which follow, we shall exhibit several explicit expressions for the vari- 
ous terms in the perturbation expansions (4), as well as comment on the accuracy of the 
approximations computed using them. 

In step two of the hybrid method, we use a subset of the perturbation coordinate functions 
{«, } determined in step one as both trial and test functions in a Galerkin type approximation. 
In particular, we seek new approximate solutions u, V, and a for it, v , and a, respectively, 
with 

N - 1 

u = Uq(x) + ^2 f>jUj(x). (9) 

J=1 

In (9), the functions {it,} are the perturbation coordinate functions determined by the 
Lindstedt-Poincare method in step one, while the {6j} represent new “amplitudes” of these 
coordinate functions. (We note that (9) satisfies conditions (lb) and (lc) for any choice of 
the {£,}.) To determine the amplitudes {£,•}, as well as the quantities u and 5, we substitute 
(9) into (la) and require that the residual is orthogonal to each iifc, 0 < k < N — 1. Also, 
we require that (2) is satisfied when u is replaced by u. Thus we obtain the conditions 

r2 7 r 

/ [? 2 u" + u + (e/a) /(a u, a v u', e)] Uk dx = 0, 0 < k < N — 1, 

Jo 

( 10 ) 

r 2 7T 

/ f(au,avu\e)u f dt = 0. 

Jo 

Equations (10) are a system of N+l equations to determine the N + 1 unknowns Su • • • ? &N-u 
V, and a. Although these equations must, in general, be solved numerically, we note that 
the solution to this system is a point in (N + l)-dimensional space, where N is reasonably 
small, while the solution to Eq (1) is a continuous function. Also, for small values of €, it is 
reasonable to assume that the unknown quantities in (10) are “close to” the corresponding 
values in the perturbation solution (4) (e.g., Sj « e 7 ). Hence, beginning with small values 
of e and then proceeding to larger values of e, good starting values are available for the 
unknown quantities in (10), which can be used with an iterative method of solution, such as 
Newton’s method. 

If we set N = 1 in (9), we find that u = Uq(x) = cos(x), while Eqs (10) reduce to 


e y2 7r ^ 

_ i / cos(x), —av sin(x), e) cos(;r)d;r, 

7r a Jo 

r 2 7T 

/ /(a cos(x), —a v sin(x), e) sin(x) dx = 0. 

Jo 


(ii) 
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Equations (11) are a system of two nonlinear equations for the two unknowns v and a. 
In particular, if we examine the solution to these equations for small values of t and let 
u — 1 + Vi t + 0(e 2 ) and a = a 0 + 0(e), we find from (11) that the equations for ?i and a 0 
are identical to Eqs (8) for V\ and a 0 - Thus, for small values of e, our hybrid solution (with 
N = 1) reproduces the perturbation approximations for u and a, at least to within terms 
which are 0(e 2 ) and 0(e), respectively. 

In the special case when / is independent of u', we see that the second equation in (11) 
is satisfied for any choice of a and hence a is an arbitrary parameter. Then the first equation 
in (11) yields the explicit expression 


v = 


1 + 


e [ 2 * „ 1 1/2 

— / f(a cos(x), e) cos (x) dx 
a Jo J 


( 12 ) 


In the next section we shall comment on the precise form of Eqs (10) for several examples, 
as well as comment on the accuracy of the hybrid approximations generated by their solution. 


4* EXAMPLES - ONE-DEGREE-OF-FREEDOM SYSTEMS 

We now illustrate some of the basic features of the hybrid solutions outlined above with 
examples of several one-degree-of-freedom systems. The small parameter e which appears 
in the nondimensional version of the oscillator Eq (la) can have several different physical 
interpretations. The examples which follow illustrate some of these different interpretations 
and hence indicate some of the classes of problems for which we feel the hybrid method 
will be especially useful. We will discuss these classes of problems, as well as some possible 
variations of the basic method, more fully in the discussion section at the end of the paper. 

Duffing Oscillator 

As our first example, we consider the Duffing oscillator which, in dimensional form, can 
be written as 

y + LO%y + ay 3 = 0, y(0) = a, y(0) = 0, (13) 

where ujq is the (linearized) natural frequency of the oscillator, a is a specified parameter, 
and the dots denote differentiation with respect to time t. We let T be the (unknown) period 
of the oscillation and then define u = 2ir/T, v = w/wo , x = tuf, and u = y/a. With these 
definitions, Eq (13) can be written in the form of Eq (la) with / = u 3 , i.e. 

v 2 u" + u + eu 3 = 0, u(0) = 1, u^O) = 0, e = aa 2 /u> q. (14) 

Thus, in this example, e can be interpreted as a measure of the amount of nonlinearity in 
the restoring force in the system. 
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For this case condition (2) is satisfied automatically. Thus, the initial amplitude a is 
arbitrary and hence can be incorporated into our expansion parameter e. The Lindstedt- 
Poincare method yields 

u 0 = cos(x), 

Ui = — [cos(x) — cos(3x)] /32, 

U 2 = [23 cos(x) — 24 cos(3x) + cos(5x)] /1024, 

( 15 ) 


21 


81 


v = 1 4- — e — 4- 

8 256 2048 


e 3 - 


6549 

262144 


37737 5 936183 6 7n 

+ — c 5 - — : e 6 + 0(e 7 ). 


2097152 


67108864 


The hybrid approximation u is obtained from Eqs (10) which, in general, will be cubically 
nonlinear in the amplitudes {£j}. In particular, for TV = 1, Eq (11) yields the solution 


1/2 


r oei 

while for N = 2, Eqs (10) yield the relations 

v = [1 +3 e/4 -3c5i/128 + 3etf?/2048] 1/2 , 

Si - e + 3e^i/4 - 9 e^ 2 /512 + 23 e<5?/16384 = 0. 


( 16 ) 


(17) 


The first equation in (17) expresses v as a function of Si and e, while the second equation 
is a (cubic) equation for as a function of e. For small values of e, we see that (16) agrees 
with (15) to within terms which are 0(e 2 ), while the solution for v from (17) agrees with 
(15) to within terms which are 0(e 3 ). 

In Figs 1 and 2, we have plotted several approximations toi/asa function of the parameter 
e. In particular, we have plotted perturbation approximations (denoted by P[A^]), where 
N, the number of terms in expansion (15), ranges from 1 to 6; hybrid approximations v 
(denoted by i/[iV]) determined from (16) and (17) with JV = 1 and N = 2, respectively; 
and numerical solutions obtained by evaluating the integral representation (3) for v. As the 
figures illustrate, the perturbation results (15) are useless for values of t above 1 (the radius 
of convergence of the perturbation series), while the one and two term hybrid results are in 
excellent agreement with the numerical results, even for values of e as large as 500. 
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Simple Pendulum 


The equation of motion of a simple pendulum can be expressed as 

9 + u >% sin (9) = 0, 0(0) = a, 9(0) = 0, (18) 

with ujI = g/L, where g is the acceleration due to gravity and L is the length of the pendulum. 
Here 9(t) is the angle the pendulum makes with the vertical. Using the same definitions of 

T, u, v, and x as in the first example and letting u = 9/a, we find that (18) can be written 
in the form of (la) as 

v 2 v!' T u + (1/e) [sin(eu) — eu] = 0, t = a. (19) 

Thus, in this case, e can be interpreted simply as the maximum amplitude of the oscillation. 
(Condition (2) is again satisfied automatically and hence the initial amplitude is arbitrary.) 
Also, since the nonlinear term in (19) can be expanded for small values of e as a power series 
in e , we see that the series in (4) involve only even powers of e. In particular, using the 

symbolic manipulation system Mathematica [15], the Lindstedt-Poincare method yields the 
results 

Uo = cos(x), 


u 2 = [cos(x) — cos(3x)] / 1 92, 

u 4 = [17 cos(x) - 20 cos(3x) + 3 cos(5x)] /61440, 


( 20 ) 


1 1 2 1 
= 1 — — e 2 + 

16 3072 


e 4 - 


23 


737280 


e 6 - 


2519 


1321205760 


e 8 + 0(e 10 ). 


This expansion converges for |e| < tt. 

The hybrid approximations u and v are obtained from Eqs (10) with (e/a) f replaced by 
[sin(eu) - eu]/e. In particular, for N = 1, Eq (11) yields 


^ = [— J sin(c cos(x)) cos(x)dx ] 1/2 

( 21 ) 

= 1 - e 2 /16 + e 4 /1536 + 0(e 6 ), as e -> 0, 

which agrees with (20) to within terms which are 0(e 4 ). In a similar manner, when the 
solutions to Eqs (10) with N = 2 are expanded for small values of e, we find that our 
expression for v agrees with (20) to within terms which are 0(e 6 ), and the solution for 6 2 
is S 2 = e 2 + e 4 /16 + 0(e 6 ). Thus, our hybrid solution u (with N = 2) agrees with the 
perturbation solution for u to within terms which are 0(e 4 ). 
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In Fig 3 we have plotted various approximations to v as a function of e. Included are 
perturbation solutions (20) (denoted by P[N)) of order 0(e 2N ) where N = 1,2,6; hybrid 
approximations (10) (denoted by H[N}) where N = 1,2 based on the first N nonzero func- 
tions {uj}\ and numerical solutions obtained by numerically evaluating Eq (3) for the case 
of the pendulum equation. The figure illustrates that the perturbation method gives good 
results for e = 9 max up to approximately e = 3 tt/ 4. For e near x, the perturbation results 
converge very slowly and give little hint that the frequency must go to zero as e -> x. In the 
region near e = 7tt/8 the H[ 2] solution seems definitely superior to the P[ 2] solution, but for 
larger amplitudes the accuracy of the H[ 2] solution is not good, though it does hint that the 
frequency is tending to zero at e = x. The H{ 3] solution (not shown) is little better than 

the H[ 2] solution. 

Large amplitude oscillations about an unstable equilibrium 

We now consider the problem of determining approximations to the frequency of “large” 
amplitude oscillations about an unstable equilibrium of a nonlinear system. As a model 
problem in this area, we consider the (non-dimensional) problem 

iP u" -u + u 3 = 0, u(0) = a, u'(0) = 0, (22) 

where u is 2 tt periodic and a is now interpreted as the nondimensional initial amplitude. 
The equilibrium state u = 0 is unstable for the system described by (22), while the states 
u = ±1 are stable. However, stable oscillations about u = 0 exist for a 2 > 2. 

To study the stable oscillations about u = 0, we consider the modified problem 

iPu" -f u + e(u 3 — 2u) = 0, u(0) = a, u'(0) = 0, 

(23) 

u(x + 2tt) = u(x), 

where e is an unspecified parameter. We now make two observations. First, we note that Eq 
(23) is of the form of Eqs (1) with f/a=u 3 -2u and hence we can construct approximate 
solutions which will be formally valid as e 0. Secondly, we note that, by setting e = 1 m 
(23), we recover our original Eq (22). Thus, our (formal) procedure will be to apply our two 
step hybrid method, as outlined above, and then set e = 1 to obtain an approximation to the 
solution to (22). In this sense, we can interpret e in this example as a type of “homotopy 
parameter. In particular, as e varies from 0 to 1, we can think of our problem (23) as varying 
from a trivial problem, corresponding to c = 0, which we can easily solve, to the problem we 
really want to solve, i.e. problem (22), corresponding to 6 — 1. 

Following the general method outlined above, we first use the Lmdstedt-Poincare method 
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on problem (23) to obtain 


uq = a cos(x), 

Ui = — (a 3 / 32) [cos(:r) — cos(3x)] , 
u 2 = (a 3 /1024) [(23 a 2 - 64) cos(;r) 

— (24 a 2 — 64) cos(3x) + a 2 cos(5x)], 


v 


1 + (3 a 2 





63 a 4 9 a 2 1 \ 3 

~25fT + l6~2J C 


6549a 8 405a 6 315a 4 15a 2 5\ 4 

262144 ~ 2048 + “512 lfT + 8/ + > 


(24) 


The hybrid approximation u is given by (9), where the amplitudes {«$_,} and V 2 are 
determined by Eqs (10) with the term (e/a) / replaced by e (u 3 — 2u). In particular, setting 
N = 1 and e = 1 in (11) we find 


£ = [3 a 2 /4 — 1] 1/2 

= (VS/2)|a| + 0(l/|a|), as |a| — + oo. 
From the integral representation (3) of u, we find for this case that 

U = [ fl2 ( 1 +sin 2 (x)) - 2]“ 1/2 dxj 

= ^1/ ' ( ! + sin 2 (a;))- 1/2 o?x] + 0(l/|a|) 

= 0.8472 |a| + 0(l/|a|), as a — > oo. 


(25) 


(26) 


Thus, from (25) we see that v = 0.866 |a| + 0(1 /\a\) as |a| — > oo, which compares well with 
the exact asymptotic result (26). 

In Figs 4 and 5 we have plotted various approximations to v as functions of the amplitude 
a. In Fig 4 it is seen from the numerical solution that the frequency goes rapidly to zero 
as a — > y/2 from above. Also in Fig 4 the perturbation expansions (24) (*P[jV] for N = 
1,2,4, 8, 16) indicate that the convergence of the perturbation solutions takes place only for 
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a very limited range of a values, perhaps \/2 < a < 1.75. The hybrid solutions (H[N\ for 
N = 1,2, 3, 4) show rapid convergence for a > 1.7 and slower convergence nearer a = \/2. 
Even so in the region just above a = y/2, H[ 4] seems much more accurate that P[ 16]. Figure 

5 shows that the hybrid solution, H [2] continues to have good accuracy up to a = 50 and 

that even H[ 1 ] is a fairly good approximation for a = 50. 

Self excited oscillations — the van der Pol oscillator 

As an example of a nonconservative system, we now consider the limit cycle of the van 
der Pol oscillator, for which f(a u,v a u', e) = (a 2 u 2 - 1) a v u' in Eq (la), i.e. 

v 2 u” + u + e(a 2 u 2 - l)vu' = 0, «(0) = 1, u'(0) = 0, (27) 


u(a: + 2 7 t) = u(x) for all x > 0. 

For this case, e is a “tuning” parameter and can be interpreted as a measure of the amount 
of a special kind of nonlinear damping in the system. 

Several terms in the perturbation expansions of u, v, and a have been reported [1,4]. In 

particular, we have 

u 0 = 2 cos(a:), 


ui = [3 sin(ar) - sin(3x)]/4, 

u 2 = [12 cos ( x ) - 18 cos(3x) + 5 cos(5x)]/96, 


v 


a 


16 3072 

1 2 1033 

2 + 96 6 552960 


e 4 + 


1019689 6 

55738368000 * 


e 6 + 


(28) 


The radius of convergence for this expansion is known to be approximately e — 1.85 (see 
Andersen and Geer [1]). 

In Fig 6 we have plotted several approximations to v, including the approximations 
P[N], obtained by including all of the terms in the perturbation expansion (28b) up to those 
which are 0(e^ +1 ), the approximations i/[./V], which are the hybrid solutions based on the 
first N perturbation coordinate functions {uj}, and the approximations of Zonnefeld [16], 
obtained using purely numerical methods. The figure illustrates that the hybrid results are 
superior to the perturbation results on which they are based, and appear to be converging 
to the numerical approximations as N increases. However, the overall quality of the hybrid 
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approximation is not as good as in the previous three examples. In the next section we shall 
discuss a general technique to improve the overall quality of the hybrid approximations and 
then apply it specifically to the van der Pol oscillator. 

5. COMBINING TWO OR MORE PERTURBATION EXPANSIONS 

In many practical applications, it may be possible to construct perturbation expansions of 
the solution u and the frequency v about more than one value of the perturbation parameter 
e. That is, it may be possible to expand both u and v as formal asymptotic series of the 
form 

u = 2>XM + 0(aS,,, + .M). 

3=0 

(29) 

v = Yh s tfie) + 0{a p Np+1 {e)), 

3 = 0 

which will be formally valid as e e p , for p = 1,2, ... ,P. Here {c*j(e)} is an appropriate 
asymptotic sequence of gauge functions and each u j can be determined completely by a 
standard perturbation method (e.g. a composite expansion of inner and outer expansions). 
For example, in the Lindstedt-Poincare method, we have constructed expansions of the form 
of (29) with c p = t\ — 0 and aj(e) = t 3 . In addition, it may also be possible to construct 
expansions of the form of (29) which will be valid as e — > e? = oo, where the gauge functions 
{a?(e)} might now include terms such as e" 9 , where q may be a positive integer or a positive 
rational number, and might also involve appropriate transcendental functions, such as log(e) 
(see e.g., [13]). 

A subset of all of the perturbation functions u* are now chosen as coordinate functions 
for the hybrid technique and an approximation u for u is sought in the form of (9), where 
each Uj is now one of the perturbation coordinate functions To determine the unknown 
amplitudes <5j, we again apply the conditions (10). 

To illustrate the application of our method when more than one perturbation expansion 
of the solution is available to us, we consider first the problem of determining the resonant 
frequency of a simple mechanical system consisting of a mass m restrained by two identical 
springs, each having natural length L and spring constant k y midway between two parallel 
walls a distance 2d apart, as discussed by Arnold and Case [3]. Thus, if we let h be the 
maximum displacement of the mass, its displacement along the centerline between the two 
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planes is described by hu(x), where 


v 2 u H + u + // u 


1 (i-(- £ 2 ”2)172] — O5 u (0) — 1? u (0) — 0, 

u(x + 2 7r) = u(x). 


(30) 


Here e = h/d, // = A/(l — A), A = L/d < 1, 1 / = 2tt/Tlo 0 , u> 0 = [2fc(l — A)/™] 1 / 2 , and 
x = where t is time. Thus e can again be interpreted as a measure of the maximum 

amplitude of the oscillation of the system. 

For small values of e, u and 1 / have expansions of the form (29), with aj = e 2 -b Using 
the Lindstedt-Poincare method, we find 


u — cos(x) + (/z/64)[cos(3a:) — cos(x)] e 2 + 0(e 4 ), 
v = 1 + (3^/8) e 2 + 0(c 4 ), as e-+0. 


( 31 ) 


For large values of e, u has a (composite) expansion in the general form of (29) with aj = e~\ 
in which the formal outer expansion, valid for |eu| > 1, must be supplemented by inner 
expansions around x — 7r/2 and x = 3tt/ 2, where u vanishes. In this case, it is easy to show 
that u is still well approximated by the first term on the right side of (31a), while 


v — 1 T fi -f- 0(c 4 ), as c — ► 00 . 


(32) 


To apply our hybrid method to this problem, we look for an approximate solution u for 

u in the form 

N 

(33) 


u = Y^Uj(x)6j, S 0 = 1, 
3 = 0 


where Uo(x) = cos(x) and the remaining (uj(x)} can be selected from either the small or 
large-e expansions of u. Substituting (33) into (10), we see that the orthogonality conditions 
become 

F j {v 2 ,6 u 6 2 ,...,6 n ) = 0, j = 0,1,..., N, 

r in ( 34 > 

Fi = h [P 2 g" + g + ^5(i ~ (1 + e » ffl) , 

Equations (34) are a system of N + 1 equations for the N + 1 unknowns V 2 and 6j, ( j = 
1,2, ... , N). For the special case when iV = 0, we can solve (34) explicitly and find 


v 2 = 1 + 


M f 2 * 

7T JO 


T - 


(1 + e 2 cos 2 (x)) 1 / 2 J 


cos 2 (x) dx. 


(35) 


If we can expand (35) for small values of e, we recover the expansion for u in (31b), while if 
we expand it for large values of e we recover the expansion (32). In Fig 7 we have plotted 
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v determined by (35) as a function of e, and have also plotted some corresponding values 
determined by the numerical evaluation of Eq (3). As the figure indicates, the agreement of 
the hybrid results with the frequencies determined from Eq (3) is excellent. 

As a second example, we again consider the van der Pol oscillator. For large values of e, 
both u and v have expansions of the form of (29) (see e.g. [13]), where the gauge functions 
now include e -1 , e~ 7 / 3 , e _3 log(e), In particular, 

1/ = ^° g(4) e -1 + 0(e“ 7/3 ), as e — ► oo. (36) 

2 7T 

The leading term in the corresponding perturbation solution for u involves several inner and 
outer expansions, which can be combined in a straightforward manner to form the leading 
term in a composite expansion of the solution. We shall denote the leading term in this 
expansion by Uoo(x,c). 

We now use u^x^t) as one of our coordinate functions in our hybrid solution (9), which 
we write in the form 

N - 1 

u — So uo “1“ ^ ^ Sj Uj + (1 ^o) (37) 

i= i 

where {uj, j = 0,1,2,...} represent the perturbation coordinate functions determined from 
the (regular) perturbation expansion about 6 = 0. Thus, our solution (37) combines N 
terms from the small-e perturbation expansion of u with the leading term in the large-6 
perturbation expansion of u. The N + 2 unknowns 6o, . . . , 6w-i, a and u are then determined 
by substituting (37) into the N + 2 conditions (10), with N replaced by N + 1 and with u n 
formally replaced by u 0 c . In Fig 8 we have plotted two perturbation approximations to v\ the 
three-term small-e approximation denoted by Po[3] and the one-term large-e approximation 
(Eq (36)) denoted by We also show three hybrid approximations: #[3,0], which is 

based on the three small-6 perturbation functions; #[0, 1], which is based on the one large-e 
perturbation function; and H[ 3, 1], which combines the information from the three small-e 
perturbation functions with the one large-e perturbation function. As the figure illustrates, 
the hybrid approximation H[ 0, 1] is a considerable improvement over P 0 c , but it gives poor 
results for e < 2. The hybrid approximation #[3, 0] is a better approximation than P[ 3], but 
it gives poor results for e > 3. Finally, #[3, 1] gives good results for very small or very large 
values of 6 and reasonable (but not accurate) results for “intermediate” (i.e. neither “small” 
nor “large”) values of e where neither perturbation expansion is accurate. 

6* SYSTEMS WITH SEVERAL DEGREES OF FREEDOM 

We now wish to study periodic solutions Q(t,e ) to systems of nonlinear equations de- 
scribing oscillating systems with several degrees of freedom, which can be expressed in the 
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form 

AQ + BQ + ef(Q,Q,Q,e) = 0. (38) 

Here e is a small, dimensionless parameter, the dots denote differentiation with respect to 
time t, Q is an M- component vector of dependent variables (where M > 1 is the number 
of degrees of freedom of the system), A and B are real, symmetric, positive definite, M by 
constant matrices, and / is an M-component vector. (Here it is convenient to require 
the components of A to be dimensionless, while the components of B and / are required to 
have the units of t~ 2 .) We assume that, for small values of e, the term ef can be expanded 
in a series of the form 

e f = € /i + I 2 + • - • (39) 

where each /, is independent of e. We also assume that the eigenvalues A j of the generalized 
eigenvalue problem associated with the linearized version of Eqs (38), i.e. 

[B-X j A]a j = 0, j = l,2,...,M, (40) 

are distinct and simple, and that the corresponding eigenvectors aj have been normalized 
so that ( Aotj , »_,) = 1 for each j. (Here ( , ) is the usual Euclidean vector inner product.) 
Consequently, the set of eigenvectors {a J5 1 < j < M } is complete and A-orthonormal. 
Then, for e = 0, Eqs (38) have periodic solutions of the form 

Q(t) = a cos(w p t - /?) a p , cu p = Ap /2 , (41) 

where a and ft axe arbitrary constants and 1 < p < M . 

Let T be the (unknown) period of a periodic solution of (38) which reduces to the solution 
(41) as e — » 0. Then we define 

UJ = v = — , x = ut - VLOpt, q(x) = Q(t). (42) 

T u> p 

In terms of these variables, Eqs (38) become 

u 2 uj 2 p Aq" + Bq + ef(q, vu p q', v 2 u\q" , e) = 0, 

q(x + 2 tt) = q(x), (43) 

where the primes denote differentiation with respect to x. 

In the first step of our hybrid technique, we again use the Lindstedt-Poincare method 
to construct perturbation solutions which are formally valid for small values of e. Thus, we 
seek solutions to Eqs (43) in the form 

q(x, e) = q 0 (x) + e q x (x) + e 2 q 2 {x) + ..., 
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qo(x) = a cos(x - f} ) a p , 
v = 1 + e V\ + e 2 i/ 2 + . . . . 


( 44 ) 


In Eqs (44), each qj is 2i r periodic in x, and each Uj is a constant. 

Substituting (44) into (43) we find that the terms which are 0(1) on the left side of (43) 
sum to zero (from the definition of <7 0 ), while the terms which are 0(e k ) with k > 1 yield the 
relations 

+ B q k = 2 v k Wp a cos(x - (3) Aa p + g k ,k = 1,2, , (45) 

where g k involves qj and i/j with j < k. In particular, 9i = ~fi{qo, w P 9o, ?o)- The 
general solution for q k will be the sum of a particular solution to Eqs (45) and an arbitrary 
multiple of qo, which is .the solution to the homogeneous version of Eqs (45). To make our 
solution unique, we impose the orthogonality condition 


r 2 7r 

/ ( Aqk , qo)dx = 0, for k > 1. 

Jo 


(46) 


We suppose at this point that k is fixed and that qj and Vj are known for j < k. If we 
now take the inner product of Eqs (45) with cos(x — /?) cf p , integrate the resulting expression 
with respect to x between 0 and 2i r, and use the periodicity condition for we find that 
the left side vanishes and hence we can express Vk as 


1 [i* 

"* = “ /„ (*’ cos < x - V dx - 


(47) 


In particular, for k = 1 in (47) we obtain 


J r2n 

1/1 = 2 7T a> 2 a Jo J l(?0 ’ Up ^ "d cos ( x ~ P) dx - 

With ^ determined by (47), we look for a solution for qk in the form 


(48) 


Jk M 

9 k = Y. cos(j(x - P)) £ Tfc,j,P «P 

i=o p=i 


(49) 


where the positive integer J k and the constants 7 ki j iP are to be determined. To determine 
these quantities, we first express the right side of Eq (45) as a trigonometric polynomial in 
the form 

Jk 

2u k u 2 p a cos (2; - /3)Aa p + g k = 9k, j cos (j{x - /?)), (50) 

. 7=0 

where the vectors {g k ,j} are independent of x and the inner product (g kt i, a p ) = 0. This 
follows from the definition of v k in (47). Substituting (49) and (50) into (45), using the 
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relations (40), and taking the inner product of the resulting expression (40) with a m , we find 
that 


A- 

J p 


(51) 


1 <m < M, (j - l) 2 + (m - p f > 0, 

while the constant 7fc,i,p is zero from condition (46). In deriving Eq (51), we have assumed 
that ^ j 2 u> 2 , unless j = 1 and m = p. 

In the second step of our hybrid method, we seek new approximations q and v to q and 
i/, respectively, where 

q - qo(x) + ^2 £j(<0<Zi( x )- ( 52 ) 

3 = 1 


In (52), the qj(x ) are the perturbation coordinate functions which appear in the expansion 
(44), while the constants {6j} represent new “amplitudes”, which must be determined. We 
note that q is r periodic for any choice of the {£_/}, since each qj is 2x periodic. To determine 
the amplitudes {^}, as well as V, we substitute (52) into (43) and demand that the residual 
is orthogonal to each of the perturbation coordinate functions qk , i.e. 


,q k ')dx = G, k = 0, — 1. 

(53) 

Equations (53) are a system of N equations for the N unknowns i, and v. In 

particular, setting N = 1 in (52) we find that q = q 0 while (53) yields the expression 


I ([? 2 wj4g" + Bq + ef{q , vu p q‘, v 2 u 2 p q" , c)] 


V = 


1 H ^ — / (/(?0) vw p q q, v 2 Lo 2 q Q ), a p j cos(z fd)dx 

■Kula Jo ' 


1/2 


(54) 


We note that, for small values of e, from (54) we can write v - 1 + evi + 0(e 2 ), where v x 
is given by (48). Thus, for small values of e, our hybrid approximation (with N = 1) agrees 
with our perturbation result up to terms which are 0(e 2 ). 


7. EXAMPLE - COMPOUND PENDULUM 


As an application of the results of the previous section, we consider two pendula, each 
of length L and each with attached mass m. They are attached in such a way as to form 
a classical compound pendulum. Letting £fi(f) and ^(^) denote the angles the two pendula 
make with the vertical, we find that 9\ and 02 satisfy the relations 


2 0i + 02 cos(0i - 9 2 ) + #2 sin (0i “ 02) + 2u>0 sin(0i) = 0 
0 2 cos(0i — 02 ) + 02 — 0? sin(0i — 02) + sin(02) = 0. 


(55) 
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In (55), u>ff = g/L, where g is the acceleration due to gravity. We now let #i — eQ\ and 
0 2 = e Q 2 to write (55) in the form of Eq (38) with e replaced by e 2 , where 


A = 


' 2 
1 



B = uj 


2 

8 



(56) 


' Q 2 [cos(e(Qi - Q 2 )) - 1] + eQl sin(e(Qi - Q 2 )) ' 
+2 e -1 to# [sin(e Qi) - e <?i] 

Qi [cos(e(<5i - Q 2 )) - 1 ]-zQ\ sin(e(Qi - Q 2 )) 

, +e -1 oSq [sin(e Q 2 ) — e Q 2 ] 


= t 


2 


■ -(1 !2)Q 2 {Qx - Q 2) 2 + QKQi - Q 2 ) - (i/3H 2 g? ' 
. — (i/2)<9x(Qi - g 2 ) 2 - g?(Gi - g 2 ) - (i/6H 2 Qi . 


+ 0(e 4 ). 


For this problem we find 


tc>l 


\J2 — V2lob, 


U) 2 


\J 2 “h 


■ 1 1 

{ 1 1 

,V2r 

“ 2 " C2 1 J 


c,=( 1/2) (2 -V2), c 2 = (1/2)(2 + V2). 


(57) 


(58) 


Thus, for small amplitude oscillations, the periodic mode corresponding to p = 1 (i.e. the 
mode with frequency u>i) represents an oscillation for which, at any instant of time, the 
pendula are always on the same side of the vertical. We shall refer to this mode as the 
“symmetric” mode of oscillation. In a similar manner, the periodic mode corresponding to 
p — 2 (with frequency u> 2 ) represents an oscillation for which the pendula are on opposite 
sides of the vertical, which we shall refer to as the “asymmetric” mode of oscillation. Then, 
setting p = 1 in our formulae of the previous section (with a — 1 and (3 = 0 in the definition of 
q 0 ), we find that we can express the perturbation solutions for the frequency 10 of oscillation 
of the symmetric mode, as well as the angular displacements 9\ and 0 2 , as 


l u 2 /wg = (2 — \/2) v 2 = 2 — \/2 + 

2 (— 102 + 71^2) 4 (470602 - 332533^) , 

£ 16 + e 14336 + 

6 (-2676206382 + 1892439443^2) 


(59) 
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0\(x) — t qi(x) = t cos(x) + 


(1 - V2) cos(x) 23(— 19 + 24^2) cos(3x) 

32 + 2688 


+ 


(233 + 1375%/2) cos(x) (812575 - 535236^2) cos(3x) 
28672 + 802816 + 


( 60 ) 


3(970497 - 751760v / 2)cos(5x) 


9748480 

0 2 (x) = eq 2 (x) = e \/ 2 cos(x) + 


] + o(£ 7 ), 


(2 — y/2) cos(x) (576 — 781\/2) cos(3x)] 


32 


2688 


+ 


[(-2750 - 233\/2) cos(x) (-1148584 + 756783\/2) cos(3x) 

I ^ . /i ! 


(61) 


28672 

(-4066960 + 3111099^2) cos(5x) 


802816 


+ 0(e 7 ). 


9748480 

It follows that the perturbation expansions for v and for the (total) energy m u>g L 2 E for 
the symmetric mode are given by 

2 (—31 + 20\/2) 4 (-61839 + 46888\/2) 

V ~ + C 32 + 6 28672 

6 (-284017367 + 197610484\/2) 


+ 1 


25690112 


+ 0(e 8 ), 


£ = 2e2 + £4 3(29 + 8V2) 

32 


, 6 (2738185 - 550288\/2) 8 

' r t 401408 + ^ 


(62) 


If 9\ (x) and 0 2 (x) are computed to order 0(e 2N *), then the energy E is conserved (in time) 
to 0{e 2N ). 

The corresponding expressions for the asymmetric mode of oscillation are the \/2 con- 
jugates of these expressions, i.e. everywhere \/2 appears in Eqs (59)-(62) it is replaced by 

-V2. 

For both the symmetric and asymmetric modes v is a monotonic decreasing function of 
E in the interval 0 < E < 6, the energy interval in which the pendulum exhibits a back- 
and-forth motion. For any given energy in this interval the frequency of the asymmetric 
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mode is higher than the frequency of the symmetric mode. Parametric plots of //(e) vs. E(e) 
as computed to orders 0(e 2 ), 0(e 4 ), 0(e 6 ), 0(e 8 ), 0(e 10 ), and 0(e 12 ), are shown in Fig 9 
for the symmetric case and in Fig 10 for the asymmetric case. These curves are labelled 
P[ 1] through P[6]. The hybrid results are useful well past the radii of convergence for the 
perturbation results. They appear to be converging to the correct values for the energy 
interval shown. 

8. DISCUSSION AND CONCLUSIONS 

Each of the examples we have considered provides some insights into the hybrid method 
which we now discuss briefly. 

The Duffing oscillator illustrates the fact that the usefulness of a perturbation expansion 
can often be limited because the expansion has a finite radius of convergence. In this case, 
the perturbation series converges only for | e j < 1, since the solution u, when regarded as a 
function of both x and e, has a singularity at t = — 1. (This follows from the fact that the 
solution to (14) is non-oscillatory for e < —1.) Thus, the perturbation expansion fails to 
converge over the range 1 < e < oc, which is most of the region of physical interest. The 
hybrid method overcomes this limitation and provides meaningful and, in this case, at least, 
very accurate approximations to the resonant frequency, even for e well beyond the radius 
of convergence of the perturbation solution (see Figs 1 and 2). 

The perturbation expansion of the frequency (and solution u) for the simple pendulum 
converges over the entire interval of interest, i.e. — n < c < 7r, and hence, in principle, could 
be used to compute approximations to the resonant frequency for any value of t in this 
range. However, the rate of convergence is very slow for values of t close to ±7 r and a great 
many terms in the perturbation series would be needed to provide a solution of acceptable 
accuracy (see Fig 3). By contrast, the hybrid approximation, based on only two terms of 
the perturbation expansion, provides an approximation to the resonant frequency which is 
essentially as accurate as the perturbation solution based on six terms for e not close to ±7r, 
and is reasonably close to the numerically computed frequencies, even when e is close to ±7r. 
Thus, we see that the hybrid method has the effect of accelerating the convergence of an 
otherwise slowly converging series. 

In the example of oscillations about an unstable equilibrium, we saw that the perturbation 
expansions depend on the parameter a (the maximum initial amplitude), as well as upon 
e. Thus, the convergence of the perturbation solutions will also depend upon the parameter 
a. In particular, for e = 1, which corresponds to the case of interest in this problem, we 
conclude from Fig 4 that the series solution for the resonant frequency converges only for a 
very limited range of values of a, this range being from a — \/2 to approximately a & 1.75. 
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This is in contrast to the hybrid solutions which appear to converge to the exact frequencies 
for all values of a > y/2. However, the rate of convergence, even for the hybrid solutions, is 
much slower for values of a close to y/2 (see Fig 4) than for larger values of a (see Fig 5). 

The van der Pol oscillator is an example of a nonconservative system, for which both the 
frequency u and the initial amplitude a of the motion appear as unknowns in the problem 
formulation. For this case, we again note (see Fig 6) that the hybrid solutions appear to be 
converging to the “exact” (numerically determined) frequencies as the number of terms in 
the approximation is increased. This convergence appears to hold even for values of e greater 
than the radius of convergence of the perturbation expansion, which is approximately 1.85. 
However, the rate of convergence of the hybrid solutions for e > 1.85 is not as dramatic 
as in the previous examples. Although we do not as yet have a good explanation for this 
behavior, we feel (intuitively) that it is due to the rather complicated form of the response 
u as e becomes large (see e.g. [13]). It simply appears that this behavior cannot be well 
approximated by the perturbation coefficient functions which have been included so far in 
our hybrid approximation. 

When perturbation expansions of the solution about more than one value of the parameter 
t are available, it has been our experience that it is usually advantagous to use at least one 
term from each of these expansions. This is illustrated in the springs-and-mass example as 
well as for the van der Pol oscillator. 

In the springs-and-mass example, the exact form of the first term in the small-e per- 
turbation expansion, i.e. Uo = cos(a:), was used in the hybrid solution, along with an 
approximation to the exact form of the first term in the large-e perturbation expansion, 
which is again cos(x). In particular, we did not include in our hybrid approximation the 
various inner expansions which should be combined with the outer expansion cos(;r) to form 
a uniformly valid first approximation to the solution when e is large. Nonetheless, the hybrid 
results are very accurate for all values of the expansion parameter e, as well as for an entire 
range of values of the parameter p (see Fig 7). This illustrates the fact that apparently one 
does not need to determine the “complete” form of the perturbation coefficient functions 
and that merely a good approximation to these functions may suffice for use in the hybrid 
approximation. 

For the van der Pol oscillator, the exact form of the large-e perturbation solution is very 
complicated, although several terms in the various inner and outer expansions necessary 
to construct this expansion have been computed [13]. However, in our approximation, we 
used only the leading term from this expansion (which is “difficult” to compute), along with 
several of the small-e perturbation coordinate functions (which are “easy” to compute) to 
form our hybrid approximation. As illustrated in Fig 8, there is a definite improvement 
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in the quality of the hybrid approximation, although the agreement with the numerically 
determined frequencies is by no means entirely satisfactory. Presumably, the accuracy of our 
approximations could be improved by including more terms from either the small- or large c 
perturbation expansions. However, we have not carried out these calculations. 

The compound pendulum is an example of the application of our method to systems 
with more than one degree of freedom. We note from Figs 9 and 10 that the perturba- 
tion expansions of the different resonant frequencies appear to have significantly different 
radii of convergence. In particular, the perturbation expansion of the higher resonant fre- 
quency, corresponding to the asymmetric mode of oscillation, has a much smaller radius of 
convergence (approximately c* = 0.092) than the expansion of the lower resonant frequency, 
corresponding to the symmetric mode of oscillation (approximately e = 1.2). Hence, the two 
perturbation expansions will be useful for calculating approximations to the corresponding 
resonant frequency over different ranges of values of energy E (approximately 0 < E < 0.17 
for the asymmetric case vs. 0 < E < 2.3 for the symmetric case). In both cases, however, 
the hybrid approximations based on only a few perturbation functions appear to converge to 
the numerically determined solutions, even for values of t well beyond the respective radius 
of convergence. We are currently developing the hybrid method in the context of general 
Hamiltonian or Lagrangian systems (and even more general systems whose solution can be 
characterized in terms of a variational principle) and are applying the methods to several 
other examples. The results of these investigations will be reported elsewhere. 

We note that the perturbation parameters e employed in the above-mentioned problems 
have a variety of interpretations, but in each case t = 0 corresponds to simple harmonic 
motion. 

While we have for reasons of space limitation only reported on frequency functions, the 
hybrid technique appears to converge to correct mode shapes as well as to correct frequencies. 

The hybrid method is more successful in some of the examples reported here than in 
others. We think it is particularly effective for large values of e in Fig 2, for large values 
of a in Fig 5, and for all values of fi in Fig 7. It seems to be less effective, as in the 
simple pendulum problem, when “directly” approaching a singularity rather than merely 
approaching a radius of convergence. 

We also note that the equations, such as (10) and (53), which determine the hybrid 
amplitudes {<$j} and v for nonlinear systems will, in general, be nonlinear and have multiple 
solutions. Some of these solutions may be ruled out on the basis that they involve complex 
numbers. Of the remaining solutions, it appears that the only ones of interest are those for 
which the 8j(e) coincide with the gauge functions of the perturbation expansion in the limit 
e _ ► 0. Therefore following a solution path starting at t = 0 seems to be an essential part of 
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the method. In our experience this leads in each case to a unique solution. 

It is interesting to note that the hybrid method provides a nice supplement to some of the 
methods and ideas presented by Van Dyke [14] for improving the usefulness of a perturbation 
expansion. In particular, Van Dyke s ideas are applicable when a fairly large number of terms 
in a perturbation expansion can be computed (usually with the aid of numerical or symbolic 
computation). The coefficients in the series are then analyzed to help uncover some of the 
analytical structure of the solution and then this information is used to recast the series 
into a different form which, in general, will be valid for a wider range of parameter values. 
In contrast, the hybrid method requires only a few terms (often only one or two terms) 
in the perturbation expansion to be computed and then uses these terms to construct an 
“improved” approximation to the solution. 

In conclusion, it appears that the hybrid perturbation-Galerkin technique is a useful way 
to enhance the usefulness of perturbation solutions to resonant frequency calculation prob- 
lems. In particular, the hybrid solutions v provide useful approximations to the resonant 
frequencies for the examples illustrated here, even for parameter values well beyond the 
radius of convergence of the perturbation solutions on which they are based. We are cur- 
rently investigating application of the method to more general Hamiltonian and Lagrangian 
systems with a finite number of degrees of freedom, as well as to systems with an infinite 
number of degrees of freedom, i.e. to partial differential equations. The initial results of these 
investigations have been very promising. 

While there are still many theoretical questions to be answered about the hybrid method, 
it is based on an intuitively plausible idea, it is relatively simple to implement, and it appears 
to provide reasonable and often very accurate approximate solutions. 
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FIG 1. Comparison of hybrid (solid lines) frequency computations for the Duffing oscillator 
with perturbation (dashed lines) and numerical (circles) computations for a small range of e 
values. 



FIG 2. Comparison of hybrid (solid lines) frequency computations for the Duffing oscillator 
with numerical (circles) computations for a large range of e values. 
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FIG 3. Comparison of hybrid (solid lines) frequency computations for the simple pendu- 
lum with perturbation (dashed lines) and numerical (circles) computations. The amplitude 
ranges from 0 to 7r. 



FIG 4. Comparison of hybrid (solid lines) frequency computations with perturbation (dashed 
lines) and numerical (dots) computations for a small range of e values. 


26 






FIG 5. Comparison of hybrid (solid lines) frequency computations with numerical (dots) 
computations for a large range of e values. 



FIG 6. Comparison of hybrid (solid lines) frequency computations for the limit cycle of the 
van der Pol oscillator with perturbation (dashed lines) and numerical (circles) computations 
for a range of e values. 
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